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We report a computer-simulation study of the equilibrium phase diagram of a three-dimensional 
system of particles with a repulsive step potential. Using free-energy calculations, we have deter- 
mined the equilibrium phase diagram of this system. At low temperatures, we observe a number of 
distinct crystal phases. However, under certain conditions the system undergoes a glass transition 
in a regime where the liquid appears thermodynamically stable. We argue that the appearance of 
this amorphous low-temperature phase can be understood by viewing this one-component system 
cis a pseudo-binary mixture. 

PACS numbers: 61.20.Gy, 61.20.Ne, 64.60.Kw 
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Many ordered materials can undergo a first-order 
phase transition that does not alter the symmetry of 
the phase. In the early work of Hemmer and Stell [1] it 
was proposed that such iso-structural phase transitions 
are to be expected if the interaction between the parti- 
cles in addition to the ordinary hard core has a "soft" 
core combined with an attractive interaction at larger 
distances. The work of ref. [H focused specifically on iso- 
structural phase transitions in materials such as Ce or Cs, 
but since then many authors have studied a wide range 
of model system that could exhibit iso-structural transi- 
tions 0, i, i, H, @, il, 8 , 1, j3, [H [H . Most of these sys- 
tems have a repulsive intermolecular potential that has a 
region of negative curvature, a feature that is known to be 
present in the interatomic potentials of some pure metal- 
lic systems, metallic mixtures, electrolytes and colloidal 
systems. The simplest example of a negative-curvature 
potential is the repulsive-step potential which consists of 
a hard core plus a finite repulsive shoulder at a larger 
radius. Systems of particles interacting through such 
pair potentials can possess a rich variety of phase tran- 
sitions and thermodynamic anomalies, including liquid- 
liquid phase transitions fsl, water- like anomalies 
and isostructural transitions in the solid region 

The repulsive step potential has the form: 
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where d is the diameter of the hard core, a is the width 
of the repulsive step, and e its height. In the low- 
temperature limit T = ksT/e << I the system reduces 
to a hard-sphere systems with hard-sphere diameter tr, 
whilst in the limit T >> 1 the system reduces to a hard- 
sphere model with a hard-sphere diameter d. For this 
reason, melting at high and low temperatures follows 



simply from the hard-sphere melting curve P = cT/cr'^, 
where c « 12 and a' is the relevant hard-sphere diameter 
(cr and d, respectively). A changeover from the low-T 
to high-T melting behavior should occur for T = 0(1). 
The precise form of the phase diagram depends on the 
ratio s = a/d. For large enough values of s one should 
expect to observe in the resulting melting curve a max- 
imum that should disappear as s — » 1 12|. The phase 
behavior in the crossover region may be very complex, as 
shown below. 

In our simulations we have used a smoothed version of 
the repulsive step potential (Eqn.[T]), which has the form: 
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+ -e(l-tanh(fco(r-a,))) (2) 



where n — 14, fco = 10. We have considered the following 
values of as', tis = 1.15, 1.35. In the remainder of this 
paper we use the dimensionless quantities: f = r/d, P = 
Pd^/e, V = V/Nd^ = 1/p. As we will only use these 
reduced variables, we omit the tildes. 

In order to get a hint about the phase diagrams of 
the system we computed energies per particle for differ- 
ent crystal ground states for ag = 1.15, 1.35. In the up- 
per part of Fig. [T]we plot a schematic phase diagram at 
T = for ijs = 1.15. It is obtained from the calculation of 
the ground state energy of the FCC (face-centered cubic) 
and BCC (body-centered cubic) crystal structures. We 
also considered several other structures for this poten- 
tial - SC (simple cubic), HCP, diamond structure, FCO 
(face-centered orthorhombic) and BCO (body-centered 
orthorhombic) but all these were unstable. 

In the lower part of Fig. [1] we show the schematic phase 
diagram at T = for as = 1.15. It is plotted from 
the calculation of the ground state energies of the FCC, 
BCC, SC, diamond and FCO structures for as ~ 1.35. 
Fig. [1] shows that there exists a range of densities where 
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FIG. 1: Ground state. The schematic phase diagrams at T : 
for Os = 1-15 and rjs = 1.35. 



the stable phase has an FCO structure whose unit cell is 
increased along the z axis in 1.6 times. This is a clear 
indication that the phase diagrams in the solid region are 
quite complex. 

To determine the phase diagram at non-zero tempera- 
ture, we performed constant-NVT MD simulations com- 
bined with free-energy calculations. In all cases, periodic 
boundary conditions were used. The number of particles 
varied between 250 and 500. No system-size dependence 
of the results was observed. The system was equilibrated 
for 5 X 10*^ MD time steps. Data were subsequently col- 
lected during 3 x l(f5t where the time step 5t — bx 10^^. 

In some region of the phase diagram for Us = 1.35 
we found that the liquid phase was thermodynamically 
stable down to the lowest temperatures that we could 
probe. In this region of the phase diagram, the dynamics 
of the system was very sluggish, forcing us to perform 
much longer MD simulations (up to 90 x 10^ MD steps). 

In order to map out the phase diagram of the sys- 
tem, we computed its Helmholtz free energy using the 
thermodynamic integration: the free energy of the liq- 
uid phase was computed via thermodynamic integration 
from the dilute gas limit [13] , and the free energy of the 
solid phase was computed by thermodynamic integration 
to an Einstein crystal [3, [3] ■ In the MC simulations of 
solid phases, data were collected during 5 x lO'^ cycles af- 
ter equilibration. To improve the statistics (and to check 
for internal consistency) the free energy of the solid was 
computed at many dozens of different state-points. All 
free-energy data were used to construct a single multi- 
nomial fit for the equation of state in every phase. The 
transition points were determined by a double-tangent 
construction (DTC). 

Fig. [2] shows the phase diagrams that we obtain from 
the free-energy calculations for two different values of 
as- Figs, [l^a-b) show the phase diagram of the sys- 
tem with as = 1.15. One can see that for the system 
with as — 1.15 there are no maxima in the melting 
curve. In a soft-sphere system described by the potential 
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FIG. 2: Phase diagram of the system of particles interacting 
through the potential (2) with iTs = 1.15, 1.35 in P — T and 
p — T planes. The asterisk in Fig. (d) corresponds to the glass 
transition. 



1 jr a face-centered cubic crystal structure has been re- 



ported 15|. However, the addition of a small repulsive 
step leads to the appearance of the FCC-BCC transition 
shown in Figs. [2Ka-b). 

Figs.[2Ic-d) show the phase diagram of the system with 
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as = 1.35 in the P — T and p — T planes. There is a clear 
maximum in the melting curve at low densities. The 
phase diagram consists in two isostructural FCC parts 
corresponding to close packing of the small and large 
spheres separated by a sequence of structural phase tran- 
sitions. 

This complexity of the phase diagram can be under- 
stood from the shape of the potential. Since the potential 
is purely repulsive, the particles tend to minimize both 
hard and soft-core overlaps. Such overlaps can be avoided 
in the low pressure part of the phase diagram. Upon in- 
creasing the density, the system crystallizes into FCC 
phase which corresponds to the close packing of the par- 
ticles at the soft core. If the pressure increases more, the 
particles begin to penetrate the soft core but still avoid 
hard-core overlaps. Hence, both parts of the potential 
are now important. The FCO phase corresponds to the 
situation when the particles are packed in X — Y planes 
in accordance with the hard core but the planes are ar- 
ranged at the soft core distance (layered structure). As 
a result a particle has 4 nearest neighbors in this struc- 
ture. If we further increase the pressure, more particles 
penetrate the soft core, but still the closed packed struc- 
tures are not favorable. Finally at very high pressures 
region we observe the FCC phase corresponding to the 
r~^^ limit. 

Interestingly, there is a region of the phase diagram 
where we have not found any stable crystal phase. In or- 
der to clarify the properties of the system in the " gap" on 
the phase diagram between the low-density FCC phase 
and FCO phase, we compute mean-squared displacement 
(MSD, Ar^(i)) and self-correlation function Fs{q,t) for 
different temperatures at p = 0.53 (Fig. [3]). As usu- 
ally observed in the proximity of liquid-glass transition 
[16|, ll7|, llSf, a bending occurs in the MSD after the 
initial ballistic regime (Fig. Elja)). A plateau appears 
at low temperatures which corresponds to the onset of 
the caging regime. At long time, the diffusion regime 
(Ar^(i) cx t) is reached, when the particles move, on av- 
erage, a distance of the order of their size. 

Another plateau is observed for self-correlation func- 
tion Fs{q,t) (Fig. [3])(b) in the time interval correspond- 
ing to the caging regime (/3-relaxation regime in mode- 
couphng theory [l^ 0, [2l|).The correlation functions 
start to decay from the plateau at times corresponding to 
the onset of the diffusive regime in the MSD (a-relaxation 
regime). In order to estimate the transition tempera- 
ture, we calculate diffusivity D as the long time limit of 
< Ar(t)2 > /6t. As predicted by Mode-Couphng The- 
ory (MCT), in the vicinity of a glass transition point 
there is the power-law temperature dependence of the 
diffusivity D cx (T - Tc)^. In accordance with MCT 
ilHIill, 7 ^ 1-75. From Fig.[3Ka) one has = 0.079 
and 7 = 2.29. 

The apparent glass-transition temperature is above the 
melting point of the low-density FCC and FCO phases 
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FIG. 3: (a) MSD and (b) self-correlation function Fs{q,t) as 
functions of time for a — 1.35 for different temperatures at 
p = 0.53. 



(see Fig. [DJd)). This suggests that the "glassy" phase 
that we observe is thermodynamically stable. This is 
rather unusual for one-component liquids. In simula- 
tions, glassy behavior is usually observed in metastable 
mixtures (see e.g. 16, [S]) where crystal nucleation 
is kinetically suppressed. One could argue that, in the 



glassy region, the present system behaves like a "pseudo- 
binary" mixture of spheres with diameters d and and 
that the freezing-point depression is analogous to that 
expected in a binary system with a eutectic point: there 
are some values of the diameter ratio such that crystalline 
structures are strongly unfavorable and the glassy phase 
is stable even for very low temperatures. The glassy be- 
havior in the reentrant liquid disappears at higher tem- 
peratures. In this density range we also found the anoma- 
lous behavior of isochores which corresponds to the neg- 
ative thermal expansion coefficient. This behavior will 
be discussed in detail in a subsequent publication. 

The behavior that we observe in this three-dimensional 
system bears some striking similarities to that observed 
in two-dimensiona.1 system with "repulsive-shoulder" po- 
tentials 22, l23|, |2J, |25|. These simulations suggest 



that such systems have a low-temperature reentrant liq- 
uid phase where particles intend to arrange in more 
or less complex, locally ordered patterns. Norizoe and 
Kawakatsu [25| find evidence for local clustering and the 
formation of percolating clusters in the reentrant liq- 
uid phase of a 3D repulsive-shoulder model. Moreover, 
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FIG. 4: (Color online). Snapshot of the disorderded phase 
for the cTs = 1.35, for p = 0.55, T = 0.1. The colors on the 
snapshots correspond to the number of the nearest neighbours 
(yellow - 0, orange - 1, dark grey - 2, red - 3, blue - 4, light 
blue - 5). 



refs. [UlilHIiil find that diffusion in this part of the 
phase diagram is very slow - for this reason, they refer to 
these percolating cluster phases as "amorphous solids". 
However, none of these studies include free-energy calcu- 
lations and hence the equilibrium phase boundaries be- 
tween solid and reentrant liquid are not known. Our free 
energy calculations support that the notion that these 
amorphous solids are thermodynamically stable. 

In the Fig. [4] we represent the snapshot of our system 
in this region for CTs — 1.35, for p — 0.55, T — 0.1. Parti- 
cles in the snapshot are colored according to the number 
of the nearest neighbours. The nearest neighbours to 
the particle X were considered to be all particles within 
the distance of the first minimum of the radial distribu- 
tion function. After reentrant melting the structure cor- 
responds to an assembly of linear "strings" of particles 
with 2 near neighbours. It is possible that the structure 
of liquid here is 3D analog of the spatial configurations 
found in [2^, this point clearly requires further study. 

In summary, we have performed the extensive com- 
puter simulations of the phase behavior of systems de- 
scribed by the soft, purely repulsive step potential ([2]) in 
three dimensions. We find a surprisingly complex phase 
behavior. We argue that the evolution of the phase di- 
agram may be qualitatively understood by considering 
this one-component system as a pseudo-binary mixture 
of large and small spheres. Interestingly, the phase di- 
agram includes two crystalline FCC domains separated 
by a sequence of the structural phase transitions and a 
reentrant liquid that becomes amorphous at low temper- 
atures. The phase behavior of systems with even wider 
repulsive steps will be discussed in separate publication. 
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